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Abstract 

This paper deals with multigrid methods for computational problems that arise in the theory 
of bifurcation and is restricted to the self adjoint case. The basic problem is to solve for arcs of 
solutions, a task that is done successfully with an arc length continuation method. Other impor- 
tant issues are, for example, detecting and locating singular points as part of the continuation 
process, switching branches at bifurcation points, etc. Multigrid methods have been applied to 
continuation problems [BK],[M]. These methods work well at regular points and at limit points, 
while they may encounter difficulties in the vicinity of bifurcation points. A new continuation 
method that is very efficient also near bifurcation points is presented here. The other issues 
mentioned above are also treated very efficiently with appropriate multigrid algorithms. For 
example, it is shown that limit points and bifurcation points can be solved for directly by a 
multigrid algorithm. Moreover, the algorithms presented here solve the corresponding problems 
in just a few work units (about 10 or less), where a work unit is the work involved in one local 
relaxation on the finest grid. 
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1 Introduction 

This paper discusses multigrid methods for computational problems that arise in the theory of 
bifurcation. Only the self adjoint case is treated here. Generally, a nonlinear equation with a 
parameter(s) is given, and the behavior of the solution as a function of the parameter(s) is required. 
This may include, for example: 

1. Continuation along solution curves (also near singular points) 

2. Detection of singularities in the marching process 

3. Locating singular points as part of the continuation 

4. Locating singular points without continuation (on fine levels) 

5. Switching branches at bifurcation points 

Continuation methods in which the problem is solved for one choice of the parameter and then 
changing it until a prescribed value is reached are commonly used. Such methods may encounter 
difficulties in case where the problem becomes singular for some choice of that parameter. Differ- 
ent continuation techniques have therefore been developed [K],[M]. In the arc-length continuation 
methods [K],[M], a new parameter, which represents an arc- length along a solution curve, together 
with an extra equation, is introduced. This makes the continuation process very robust at regular 
and at limit points (defined later). 

Detection of singularities as part of the marching process may be required. A common way 
of doing that is by monitoring the sign of the determinant of the linearized discrete operator [K]. 
Such an approach works well if the discrete problem has a singularity whenever the continuous 
one has. However, in non trivial cases this may not be so. Bifurcation on the continuous level 
may become an imperfect bifurcation on the discrete level; that is, the different branches which 
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intersect on the continuous level do not intersect on the discrete level, while they approach each 
other, approximating better and better the continuous branches, as the discretization parameter 
goes to zero. In such cases the determinant may not change sign, and as a result some bifurcation 
points (of the continuous level) may be skipped. In this paper detection of singular points is done by 
monitoring (on the coarsest level) the eigenvalue which is the closest to zero. If a small extremum is 
detected in the behavior of the eigenvalue, further checks are done in order to ascertain the existence 
of a singular point of the continuous level in that vicinity. In case the extremum approaches zero, 
as the discretization parameter tends to zero, a singularity has been detected. The advantage of 
this method is that it can detect every singularity of the continuous system. 

Once a singularity is detected, its exact location may be required. Several approaches for 
locating singularities exist. These methods are basically root finders for the determinant of the 
linearized system and involve iterations in which at every step the nonlinear problem has to be 
solved. Such methods require 0(n 3 ) operations for each step (in calculating determinants). The 
method suggested here uses a set of equations that defines the singular point and solves these 
equations directly, without involving iterations within iterations as other methods do. This can 
lead to algorithms that need O(n) operations in order to locate singular points. The importance of 
this becomes clear when curves (in the parameter space) of singular points ( envelope of singular 
points) are required. Such envelopes of singular points can be defined as solutions of a new system 
of equations. A continuation method can then be applied to this new system, resulting in huge 
savings in the computational work. 

Having the basic tools for continuation and other objectives, a fast way of solving the resulting 
system is desired. This can be done very efficiently using multigrid methods. In [BK] a multigrid 
continuation method that uses arc-length continuation is described. It is very efficient at regular 
and limit points, but may fail near bifurcation points. In general, the difficulty with standard 
multigrid continuation methods near singular points arises because of the eigenfunctions corre- 
sponding to nearly singular eigenvalues. The coarse levels approximate these components very 
poorly, resulting in a slow convergence rate or even divergence of the multigrid cycle. In [BT] a 
method for treating nearly singular problems has been proposed. That method treats separately 
the poorly approximated components, which turn out to be the nearly singular components. It 
begins by identifying them and then introduces a correction to the usual coarse grid equations in 


4 


such a way that makes them approximate all components well enough for accelerating finer grid 
convergence. The amount of work involved in computing the nearly singular components to the 
desired accuracy is very small, making the cost of the modified algorithm essentially that of the 
standard one, except in linear problems which are extremely close to singularity. In that case the 
accuracy of nearly singular components is improved as the algorithm proceeds in getting better 
and better approximations for the original function to be found. 

The paper is divided into two parts. The first, which includes sections 2-4 , discusses basic ideas 
concerning the treatment of bifurcation problems. These are arc length continuation, the detection 
and location of singularities (or envelope of singular points), and switching branches at bifurcation 
points. This part of the paper except part of section 3.2, which is believed to be original, is given 
for completeness of the presentation of the subject . 

The second part is devoted to the description of various multigrid methods. It begins in section 
5 where standard multigrid methods are described. This includes the multigrid method for contin- 
uation problems at regular and limit points. This method is basically the frozen tau method [B]. 
This section describes also a method for locating limit points directly by an appropriate multigrid 
algorithm. 

In Section 6 a description of advanced multigrid methods designed to work efficiently near 
singularities is given. These methods involve modifications of the coarse grid equations and are ex- 
tensions of ideas developed in [BT]. Locating bifurcation points directly by an appropriate multigrid 
method is also decribed there. 

In section 7 the effectiveness of some of the algorithms presented in the paper is demonstrated. 
It is shown that limit points are located using the 1-FMG algorithm for an appropriate system 
of equations defining these points. Similar results using the 2-FMG algorithm are obtained when 
bifurcation points are to be located. Also the effectiveness of the new continuation methods near 
singular points is demonstrated and compared with a standard continuation multigrid method. 

2 Arc Length Continuation 

Consider a nonlinear problem, say, in the form 


L{u, A) = 0 


( 2 . 1 ) 
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where L : X x Z —> M, Z is the set of real numbers and M is some Hilbert space and such that 
L u (u, A) is a self adjoint operator. Smooth branches of solutions 

{r a & : [u(s), A(s)], s a < 8 < *i} (2-2) 

where u(s) € Z,X (s) 6 Z are required. Here the parameter s is arbitrary, and we refer to it later 
on as the arc-length parameter. 

The standard approach of using A as the continuation parameter is basically the following. 
Given a solution uq at Aq one seeks a solution at Ao + SX by using Newton’s iteration 


lW8uW{ Ao + 5A) = -L^ n > 1 (2.3) 

where: 

4 n) = M« (n) (Ao + «A), Ao + 6 A) (2.4) 

L (n) = L(uW( A 0 + $A), A 0 + SX) (2.5) 

= uo (2.6) 

« (n+1) (Ao + SX) = uW(A 0 + £A) + £u( n >(A 0 + SX), n > 1. (2.7) 


This continuation method encounters difficulties when L u becomes singular . A method that 
circumvents such difficulties to some extent is the arc length continuation. In this method (2.1) is 
replaced by 

L(u(s),A(s)) = 0 (2.8) 

i\T(u(s), A(s),s) = 0 (2.9) 

where NtZxZxZ^Z and s G Z is the independent parameter on the arc of solutions [K],[M]. 
Several choices for N are possible. Throughout this paper N is defined by 

N(u(s),X(s),s) = 0<u- uo,«o — «i > +(1 - 0)(A — A 0 )(A 0 - Ai) - 

(8 - S 0 )^||«1 - «o|| 2 + (1 - 0)1 Ai - Ao I 2 (2.10) 

where so and are two previous points on the branch of solutions, u = u(s),«o = «(so),«i = 
u(si), Ai = A(«i), 0 < 0 < 1 and < ., . > denotes the inner product in M. This equation is an 
approximation to a choice given in [K] and does not require the computation of the tangent to the 
solution curve. 
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Given (u(s 0 ), A(s 0 )) and (ti(si),A(si)), a process of solving for (u(s),A(s)) proceeds as follows: 

..0/_\ _ , ( g ~ a l) \ m foil! 


u°(s) = u(s 0 ) + ( Sl - 7j HSl) ~ 

A°(s) = A(s 0 ) + ^ _^) t A (*i) - A (s 0 )] 
u(")(s) = «( n - x )(s) + 5«( n - x ) n > 1 
AW( S ) = A( n - X )(s) + SA( n - x ) n > 1 
where (<5u( n ),6A( n )) is the solution of the equation 

( 4 n) \ ( Sul”) \ ( -2» ' 


N { u n) N[ n) 


-2>) 

~N( n ) 


( 2 . 11 ) 

( 2 . 12 ) 

(2.13) 

(2.14) 


(2.15) 


and iffl, lS"\ lS n \ Nu ^ , N[ n \ are defined in a similar way as in (2. 4), (2.5). A multigrid 

method for solving equations (2.8),(2.9) at regular and at limit points is described in section (5.2). 
An extension of that multigrid algorithm for solving near bifurcation points is given in section 6. 

The role of 9, in the definition of N , is to enable inexpensive marching through limit points 
where the curve has a high curvature. Note that the geometric meaning of (2.10) is that the angle 
between two successive changes in the solution, in the (u, A) plane, is less than tt/ 2 ( assuming 
a > $o). In case the step size is too large, for example, in encountering an angular limit point, 
this may not be the case, and the augmented system (2.8) (2.9) may not have a solution unless s 
is close enough to so- With a proper choice of 9 (usually either 9 — 1 or 9 = 0) marching through 
such points is possible, without decreasing the step size . 

A Continuation Process 

For completeness of the presentation, a full continuation process is described here. Let 9 = 0.5. 
Step 1: 

• Set so = 0, A = Ao,uo = 0. 

• Solve L(u o, A) = 0, for tio , keeping A fixed. 


Step 2: 

• Set A = Aj, tii = uo . 

• Solve L(u i, A) = 0, for ui , keeping A fixed. 

• Compute si from the equation 

(si — so) 2 = 0||tti — «o|| 2 + (1 - 0)| Ai - A 0 | 2 . 
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Further Continuation Steps 

• Choose As (for example si — so ), Set s = si + As and 

“ = “i + ££[<*1 - «»l 
a = AI + ££[A, - A„) 

• Solve the following equation for (u, A) 

L(u, A) = 0 

N(u,X,s) = 0 

• Set so = $i, ®i = s, «o = tti, ill = u 

3 Singularities 

3.1 Detection of singularities 

A singular point is a point (u*,A*) on a branch of solutions for which L u (u* , A*) is singular. 
A common way for detecting singularities is to check for a sign change in the determinant of the 
linearized system [K] . This, however, cannot detect bifurcation from eigenvalues of even multiplicity. 
Also it cannot detect bifurcation of the differential equation which became imperfect bifurcation 
on the discrete level, since the determinant may not change sign along discrete solution curves in 
such cases. In fact singularities may disappear upon discretization. 

In this paper a different approach is taken. Let (u*,A») be a singular point, and let <p be a 
corresponding singular eigenfunction, i.e., 


L u (u*, A*)y? = 0. 


In the vicinity of the point (u*, A») 

H 2 {s) = min 
(P# o 


< L u (u(s),A(s))»?,y? > j 2 
L(u(s),A(s)) = 0, 


subject to 


(3.1) 


(3.2) 

(3.3) 


is close to zero, and it vanishes at (u, A) = (u., A„). Detecting singular points is done by computing 
approximately p? (s) on the coarsest level. Once an extremum has been found, and it is small 
enough, further checks are needed in order to decide whether a zero for n(s) exists on the continuous 
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level. For the continuous level the extremum of ^i(s) 2 is zero at singular points, while it may not be 
so on the discrete level. However, it approaches zero as the discretization parameter goes to zero. 

In detecting singular points we distinguish between limit points, bifurcation points, and imper- 
fect bifurcation points. In the last case L u may not become singular but almost singular; that is, it 
has a very small eigenvalue. This case is very important since, as mentioned above, a bifurcation 
in a differential equation may become an imperfect bifurcation for the discrete system, and this is 
very likely to happen (See section 7) . 


Limit Points A point (uq, Ao) is called a limit point if 


d«mA/(L°) = 1 (3.4) 

codimR(Lu) = 1 (3.5) 

Ll * Z(L° U ), (3.6) 

where £(L°) are the null space and the range of = L u (u o,Ao) respectively. These 

conditions imply that 

Ao = 0. (3.7) 


Bifurcation Points A bifurcation point is a point where two branches of solutions have a non 
tangential intersection. At such points 


d»m^/(L°) = codimZ(L °) = m 

Ll e *(L°). 


(3.8) 

(3.9) 


Imperfect Bifurcation Imperfect bifurcation is a situation in which a problem with a bifurca- 
tion has been perturbed and as a result the different branches do not intersect any more, but they 
come very close to one another. Two cases are possible: 

(i) The different branches consist of regular points only (in that vicinity ) 

(ii) At least one of the branches has a limit point. 

The next section describes how locating singularities is done and how a distinction is made once 
a singularity has been found. In this paper only singularities of the continuous level are of interest. 

Limit point and bifurcation points (perfect and imperfect) are showm schematically in figures 


1 and 2. 
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3.2 Locating Singularities 


Having detected a possibility of a singular point using the coarsest level, further checks (to confirm 
the existence of a singularity) and its exact location may be required. In the vicinity where /i(s ) 2 
has a minimum, the following minimization problem may be used to accurately locate a singularity 


min m 2 subject to 
L(u,X ) = 0 
L u {u,X)<p = n<p 

IMI 2 = i. 


(3.10) 

(3.11) 

(3.12) 

(3.13) 


This minimization problem can serve for locating either limit points or bifurcation points. Together 
with /x 2 (s) also < L\ } <p > is needed in order to decide about the nature of the singularity located, 
i.e., whether it is a limit point, a bifurcation point, or an imperfect (discrete) bifurcation. At 
bifurcation points < L\,<p >= 0, and at limit points < L\,<p > ^ 0, while at both /z(s) = 0. When 
a discretization is used to approximate a singular point, the two quantities /i* and < L^, <p h > are 
available. The nature of a singularity is determined as follows: 

IF fJL* k /► 0, as A — ► 0 THEN no singularity exists 

ELSE 

IF < L\ y (p h >— ► 0 as h — ► 0, a possible bifurcation point has been detected 

ELSE a limit point has been detected. 

The location of the singular point, if detected, is obtained from the limit of the sequence (tij , Aj) 
which minimizes the discrete analogue of (3.10)-(3.13). Note that the decision about the existence 
of a singular point can be made in general only by solving a sequence of discrete problems, since a 
perfect bifurcation may become an imperfect one upon discretization. 

Sometimes it is known in advance that only limit points are involved. In such cases they can 
be located using the following system of equations, once the vicinity of the limit point is reached. 


L(u, A) = 0 
L u (u,\)<p = 0 

IMI 2 = i- 


(3.14) 

(3.15) 

(3.16) 
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In a vicinity of a limit point the above system is regular. To see this differentiate the equation 
L(u, A) = 0 twice with respect to the arclength parameter s, getting 

LJti + [L uu iiu + 2L U \ uA + L»AA] + L \ A = 0. (3.17) 

At a limit point L\ £ £(L U ),A = 0, hence it follows that L Uil uu ^ 0. Since u = a (p for some 
real a, it follows that L uu <p<p ^ 0. This implies that the linearized equations of (3.14)-(3.16) are 
not singular. A Newton iteration therefore may be used for solving these equations. A different 
method, that combines multigrid ideas, for solving the above system is described in section 5.3. 

3.3 Locating Envelopes of Singular Points 

When nonlinear equations are given in terms of more than one parameter, say for example, A and 
p, a curve of singular points (for example limit points), in the (A,p) plane may be required. A set 
of equations that define such curves is given by 


min p? subject to 

Kp.*ip) 

(3.18) 

L(u,X,p) = 0 

(3.19) 

L u (u,X,p)<p = p<p 

(3.20) 

II 

h- 1 

(3.21) 


These equations are solved once the vicinity of a singularity is reached. One possibility in solving 
these equations is to find for each p a A for which a singularity exists. This however may fail 
near limit points in the (A, p) plane. To circumvent this difficulty an equation like the arc length 
equation is added, but this time in the (A,p) plane. That is, an envelope of singular points is found 
by solving (3.18)- (3.21) together with 


N(X,p, s ) = 0(A-Ai)(A! -A 2 ) + (1-0)(?-pi)(pi - pi) - 

( S - 8 0 )y/9(Xi - A a )> + (1 - 0)(P! - p 2 )2 = 0, (3.22) 


where s is the arclength parameter, and (Ai, pi), (A 2 , P 2 ) are two successive points on the curve 
that have been computed previously. The first two points on such a curve are found, for example, 
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by fixing p and trying to find a singular point in A. This is done using methods described in 
the beginning of the section. Other points are solved for using (3.18)-(3.22). The nature of the 
singularity is found by computing < L\ ,<f>>. A simplification occurs when only limit points are 
involved; then the equations are 


L(u,\,p) = 0 

(3.23) 

L u (u,\,p)<p = 0 

(3.24) 

N(X,p,s) = 0 

(3.25) 

IMI 2 = i. 

(3.26) 


4 Switching Branches 

Once a possibility of having a singular point has been detected, a search for new solution branches 
may be required. Several methods for doing this are possible [K], Two of these are described 
here. The first is useful only when the discrete equations have a bifurcation point whenever the 
continuous problem does. The second can be applied in more general cases. 

4.1 Method I (perfect bifurcation) 

In this method switching branches at bifurcation points starts by finding directions for the new 
branches. Given the singular eigenfunctions {<pi , . . . , <p m } at a bifurcation point (u* , A*), the solu- 
tion for any of the branches going through the bifurcation point is approximately of the form 

m 

( u * + E a i t P)> A* + ar°) (4.1) 

/= o 

in the vicinity of that point, where <po is the solution of the system 

L° u <p 0 + L° x = 0, <<pj,<p o >=0 1 <j< m. (4.2) 

The superscript ® on operators indicates here that they are evaluated at the bifurcation point. The 
m + 1 scalars {a 0 , . . . , a m ) are a solution of the quadratic system of equations [K] 

m m m 

EE aijic&j&k + 2E bijOtjato + Cj-Oq = 0, 1 < t < m; 

;=!*=! j = 1 


(4.3) 
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where 


m 



j=0 


a ijk =< L° uu (pj(pk, <pi >, 1 < i,j, k < m] 

b i3 - =< {L° uu <po + Ll x ]<p,-, <pi >, 1 < i,j < m; 

Ci =< [L° uu <pq<pq + 2L° uX (p 0 + L° xx ],<pi >, 1 < t < m; 


(4.4) 


(4.5) 

(4.6) 

(4.7) 


and e is sufficiently small. 

Having a solution to (4.3)-(4.4), a new point (u,A) on the new branch can be found by solving 
the system 


L(u, A) = 0 (4.8) 

TO 

YL l< u - u o ,<Pj >| 2 = e 2 (4.9) 

i - o 

starting with the initial approximation 


m 

u = u* + J2 a i ( Pj (4.10) 

3=0 


A — A* + ao- 


(4.11) 


When this has been done a continuation that uses arclength for marching on the new branch of 
solution can be started. 


4.2 Method II (perfect or imperfect bifurcation) 

In this method solutions are required on some plane "parallel” to the tangent but displaced from the 
bifurcation point (either a perfect or an imperfect) in a direction "normal” to the tangent. Let (ip, 6) 
be a normal to the tangent (u(s), A(s)). Starting from the initial approximation (uo + tip, Ao + eS), 
we look for a solution of the system 


L(u, A) = 0 


(4.12) 


< u - u 0 ,ip> +(A - A 0 )$ = e[||V>|| 2 + (S) 2 ], 


(4.13) 


where e is chosen to be not too small in order to avoid going back to the previous branch* In fact, 
€ can be estimated from information obtained during the computation for locating the singular 
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point. Let (u* , Aj) be the solution of the discrete version of (3.10)-(3.13). Let a and c be constants 
such that — u 2 **|| 2 + | A* — A 2fc | 2 < ch a . Then e = 4 cjlrx should be large enough (it is about 

4 times the distance from (uj, Aj) to the differential singular point). For simple bifurcation cases tp 
can be approximated from differences of uj on successive levels. The multiple bifurcation case will 
require the computation of the nearly singular subspace in order to define more than one direction 
for a new branch to bifurcate. 

5 Standard Multigrid Methods 

In this and the next section we describe multigrid methods for the problems mentioned in the previ- 
ous sections. These include path following at regular and singular points, detection of singularities 
as part of the continuation procedure, locating accurately singular points, and switching branches 
at such points. Continuation procedures in the formulation of this paper require the solutions for 
two different values of A’s . These are obtained using a basic multigrid method. 

5.1 Basic Multigrid Method 

Local relaxation and a coarse grid correction form the basis of a multigrid cycle. 

5.1.1 Local Relaxation 

The main element of a multigrid method is a relaxation method, usually a local one, which has good 
smoothing properties. For linear scalar problems Gauss-Seidel is the first candidate (for systems 
a Distributed Gauss-Siedel or Collective-Gauss-Seidel may be the analogue). When non-linear 
equations are involved, a Newton-Gauss-Seidel relaxation is used. In this relaxation the equations 
are scanned in some order where for each equation one variable is changed using one Newton step 
for the nonlinear equation that this variable satisfies (keeping all the rest fixed). The use of only 
one Newton iteration is usually enough to guarantee good smoothing. 

Gauss-Seidel relaxation does not converge for indefinite problems, though it may have very 
good smoothing properties. It should be used, therefore, on fine levels, as long as the divergence 
of the smooth components is slow enough . On coarse levels a different relaxation is needed. 
Kaczmarz relaxation for linear problems and a Newton-Kaczmarz for nonlinear problems are used 
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in the indefinite case on coarse levels. In the Kaczmarz relaxation each variable gets a correction 
proportional to its coefficient in the equation such that the equation is satisfied. In a Newton- 
Kaczmarz relaxation each variable gets a correction proportional to its coefficient in the linearized 
equation such that this (linearized equation) holds. 

The problems considered in this paper are nonlinear; therefore, the multigrid method described 
is the Full Approximation Scheme which is suited for such problems. Here the coarse grids approx- 
imate the full solution rather than corrections as in the Correction Scheme (CS). See [B] for more 
details. 


5.1.2 Multigrid Cycle. Full Approximation Scheme (FAS) 

Consider a sequence of grids Clk(k < M) with mesh sizes h * satisfying . Suppose on 

each grid an operator L k is given in such a way that L k {k < M) is an approximation to L * +1 and 
the fi* grid equations are 

L k U k = F k . (5.1) 

Assume also that interpolation operators /*_!, from coarse to fine grids , and restriction operators 
I k ~ 1 9 j£~ x , from fine to coarse grids, are given. 

Given an approximate solution u k to equation (5.1), the multigrid cycle for improving it is 
denoted by 

u k «- MG(k,u k ,F k ) (5.2) 

and is defined recursively as follows: 

IF k=l THEN solve (5.1) by enough relaxations (to achieve a desired accuracy) 

ELSE 

• Perform v\ relaxation sweeps on (5.1), starting with u k and resulting in a new 
approximation ti*. 

• Starting with u k ~ x = I* -1 *** make 7 successive cycles of the type 

it*' 1 4 - MO(k - l,t i*-\ l£~ l (F k - L k u k ) + L fc- 1 7£ - 1 u fc ). 

• Calculate ti* = u k + — I^~ x u k ]. 
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• Perform i >2 relaxation sweeps on (5.1) starting with u k and yielding u k , the 
final result of (5.2). 


The cycle with 7 = 1 is called V cycle or V ( 17 , ^ 2 ) and the one with 7 = 2 is called W cycle or 
W{u u u 2 ). 

5.1.3 Full MultiGrid (FMG) Algorithm 

In order to obtain full efficiency, the first approximation on a given level is obtained from a solution 
of the same problem on the next coarser level, which itself has been calculated in a similar way. 
The resulting algorithm is called Full MultiGrid (FMG) and is described next. 

Let nJ.j be an interpolation operator (usually of higher order than I k _i)- Given the problem 
L m U m = F m , the N-FMG solution of that problem is: 

Initial setup: 

Set u M = 0. 

FOR k = M- 1, ,1 DO: 

•F k = I k +i F k+1 
•u k = 0 . 

N-FMG Algorithm 

Calculate u 1 the solution of (5.1) for k=l by several relaxations. 

FOR k = 2,...,M DO: 

• Calculate u* <— u k + Il*_ 1 u fc_1 . 

• Perform the cycle u k <— MG(k,u k , F k ) N times. 


5.2 MultiGrid Continuation Method (regular points) 

In this section the multigrid solution of the following equations is described 

L(u,X) = F (5.3) 


JV(u,A,s) = /, 


(5.4) 
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where N is given by (2.10). 

Having solved L(u, A) = F for two different values of A, A 0 and Ai, using the method of section 
(5.1), the actual continuation method, also using equation (5.4) can be done. A relaxation scheme 
that will also treat equation (5.4) is needed. 

5.2.1 Relaxation for Continuation Problems (regular points) 

The relaxation of a system like (5.3)-(5.4), which includes global parameters as unknowns, is done 
in two steps. The first step is the relaxation of equation (5.3) for u only keeping A fixed, designed 
to have good smoothing properties. The second step changes A and u. This will be referred to as 
the global step since it involves changing the global quantity A and a global component of u, namely 
4> = ui — uo- Away from singularities the local relaxation process converges reasonably well on 
coarse grids for all components involved. The global step in this case is designed such that equation 
(5.4) holds. This could be done by changing A only; however, this might have a bad effect on the 
other equation (5.3). A change that satisfies (5.4) and also takes into account equation (5.3) must 
allow for the change of a quantity other than A. This quantity is the (f> component of u. Using the 
above <f> the global step is given by 


u + /3<f> (5.5) 

A < — A + 8 (5*®) 

where /?, S is the solution for the following equations 

< L{xi + (3<f> , A + 5), <f> >=< P, <f> > (5,7) 

N(u + p<j>,\ + 8,s) = f. (5.8) 


Note that the global step is designed such that the residuals of (5.3) are zero in the <f> direction. 
The solution of the above system is approximated by the solution of the linearized equations, i.e., 

P < L u (u,\)<f>,<f> > +8 < L\(u,X) y <f> >=< F - L(u, A),<£ > (5.9) 


p < (f> y ui - u 0 > +6 ( Ai - A 0 ) = / - A(u,A,s). 


(5.10) 
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Furthermore, < L u (u,X)<f>,<f> > and < L\{u,X),<j> > are approximated using differences, 

< L u (u,X)<f>,</> >« j < L(u + e<f>, A) - L(u, A), <f> > (5.11) 

< Lx[u,X),<j> >« - < L(u, A + e) — L{u,X),<f> > . (512) 

Both approximation steps are justified since regular points are considered here. 

Let ki be the finest level on which the global step is used; (usually fcj = 1 or 2), the relaxation 
of (5.3)-(5.4) is then 

IF k < ki THEN Perform local relaxation for u in equation (5.3) 

ELSE Perform a local relaxation for u in equation (5.3) followed by a global step. 

5.2.2 Standard Full-MultiGrid-Continuation (FMGC) Method 
Given a sequence of problems 

L k {u k ,X) = F k (5.13) 

N k (u k ,X,s) = G k (5.14) 

defined on a sequence of grids fl* with mesh sizes h*.(fc < M) and such that coarse grid operators 
are approximations to finer ones as before. We describe now the continuation multigrid method. 

Step 1: 

• Set so = 0, A = Ao,u m = 0. 

• FOR k = l,...,M -1 DO: 

F k = l£ +1 F fc+1 

u k = 0 . 

• Perform N-FMG algorithm for equations (5.11) with k = M, having A fixed. 

Call the resulting approximation on level M tig*. 

Step 2: 

• Set A = Aj, u M = Uq 1 . 

• FOR k = M — 1,. .. ,1 DO: 

t ,fc _ t* „*+i 

u — i fc+l U 0 

F k = I k +1 F k+1 + L k (u k , A) - I k +1 L k+1 (u k+1 ,X) 
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• Perform N-FMG algorithm for equations (5.11), having A fixed . 

Call the resulting approximation . 

• FOR k=M-l, . . . , 1 DO: 

..k _ Jk k+ 1 
U 1 ~ i *+l u l 

F k = l£ +1 F k+1 + L k (u k ,X) - I k +1 L k+1 (u k+1 ,X) 

G k = 0 

• Compute si from the equation 

JV 1 (uJ, Ai,si) = 0. where N 1 is defined by (2.10) replacing k by 1. 


N-FMG C Algorithm 

• Choose As (for example si — so ), Set s = Si + As 

• Perform N-FMG cycle for (u k ,X k ) in equations (5.13)-(5.14) having s fixed, where the 

relaxation used is now the one of sec (5.2.1). 

Call the resulting approximation on level M (u M , A). 

• FOR k=M-l,...,l DO: 

u k = I£ +1 u k+1 

F k = I k +1 F k+1 + L k (u k ,X) - I k +1 L k+1 (u k+1 ,X) 

G k = l£ +1 G k+1 + N k (u k ,X,s) - N k+1 {u k+1 ,X,s) 

• Set sq Si — s. 

• FOR k=l, ...,M DO: 
u o = u i> u \ = u>c - 

Observe that the right hand side of the coarse grid equation when solving first for such levels is 
not just the restriction of the finest level right hand side. There is an additional term. This term 
is called in [B] the r term, and it represents the relative truncation error of the two grids involved, 
a fine and a coarse grid. The method of using this term from a previous step in a continuation 
process is called frozen-tau method. It has been introduced in [B] and has been used in [BK], The 
main idea is as follows. After the solution has converged for some A, the right hand sides of coarse 
grids are no more the restriction of the right hand side of the finest grid. Also since the quantities 
change continuously, the additional term on coarse levels changes continuously. If r of the new 
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point on the solution branch was given on the coarse level, the accuracy of the fine grid solution, 
on the coarse grid, could be obtained. By using the previous r term a good approximation to the 
final one is used , yielding a very efficient process. 

Another way of looking at introducing this term is by attempting to solve at the new A for the 
increment in the solution. Using FAS this results in the same equation as we have. 

An algorithm that is very similar to the N-FMGC has been described in [BK], The main 
difference between that algorithm and the one here is that the arc-length equations are assumed 
given on the finest grid in our case while it was used only on the coarsest level in [BK]. Another 
difference involves the method of solving the coarsest grid equations. In [BK] a Newton iteration 
for the augmented system has been used. 

5.3 Multigrid Method for Locating Limit Points 

Once a limit point has been detected by observing the behavior of A (s), for example, it may be 
required to actually locate it. This is done by applying an FMG algorithm to a discrete version of 

(3.14) -(3.16) once the continuation process has reached the neighborhood of the limit point. The 

relaxation process is basically what needs to be described. The rest is a standard FMG method for 
equations (3.14)-(3.16) using the FAS formulation, where U k in (5.1) stands here for (it*, A). 

Relaxation 

A local process is used to smooth the error in both u and <f> by relaxing the discrete version of 

(3.14) for u and (3.15) for <f> (keeping (u, A) fixed). This local step is employed on all levels. On the 
coarsest level also a global step is employed. It begins by updating the norm of <j> by multiplying 
it by a proper constant, to satisfy the norm condition on this level, followed by 

u<-u +04, A A + 5 (5.15) 

such that (/ 3 , 8 ) satisfies 

ft < L u (u, <j> > + 8 < >=< F — L{u, A), <f> > (5.16) 

ft < L uu (u,\)<f>tf>,<f> > +8 < L u \<i>,<j> >=< G - L u {u,\)4>,4> > . (5-17) 

Here F and G are the right hand sides for the corresponding coarsest grid equations, respectively. 
These equations are obtained as in section (5.2.1) upon linearization in (ft, 8). Such a linearization 
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is justified since the system (3.14)-(3.16) is not singular in a vicinity of a limit point as shown in 
section 3.2. Note that since the operator L u is singular at a limit point a local relaxation will be 
slowly converging for the 4> component in u. This, however, is taken care of by the global step 
which changes u exactly in the direction of slow convergence. 

When the envelope of limit points is required, the above procedure is combined into an FMGC 
algorithm, where here the global step has to update also the parameter p. 

The effectiveness of the multigrid algorithm for locating limit points is demonstrated in section 
7.1, where limit point for the Bratu problem is computed to the level of discretization errors using 
1-FMG algorithm. 

6 Advanced Multigrid Methods 

6.1 Difficulties Near Singularities 

The methods described in the previous section are very efficient away from singular points. Note 
that a limit point is not a singular point for the augmented system of equations; hence no diffi- 
culty is encountered there. In the neighborhood of singular points of the inflated system, such as 
bifurcation points, difficulties with the standard methods may occur. In [BT] a detailed discussion 
of the difficulty together with an algorithm for overcoming it is presented. For completeness the 
explanation given there is stated here. 

The first difficulty, which is the less serious of the two that are discussed, is the slowness of 
local relaxation processes near singularities. The global step as mentioned in section 5.2.1 does not 
always improve the rate of convergence. In order for this to happen the function <j> used in the 
global step should be an approximation to a slowly convergent component in the local process; in 
case few such components exist, more than one function should be used. These components turn 
out to be the almost singular eigenfunctions of the operator L u (tz, A). The global step should use 
as many functions as there are components that are slow to converge. It is shown below that the 
closer one is to singularity the more accurate the eigenfunctions have to be approximated. 

Consider a linear equation Ax = 6, where A is a symmetric operator, and let (^y, Ay) be 
its eigenfunctions and eigenvalues, respectively. Let <f > i be slowly converging in a local relax- 
ation process. Assume for simplicity that the residuals consist of that eigenfunction only, i.e., 
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b — Ax = . Let the global step be x «— x + aip where ip = <p 1 + ay^y and suc h that 

a minimizes ||A(* + aip) — i|| 2 (in the general case) or 

a minimizes < A{x + aip),x + aip > —2 < b,x + aip > in the case A is positive definite. 


These global steps imply a = • < <x^*AV’> > an< * Q = < <a1(>*$>' respectively. For our special 
residuals it means a = , , - v * w — or a = , , v>n x ,,. > , respectively. In both cases it is 

seen that a good convergence implies that ay, (1 < j < m) should be small, and the closer Ai is to 
zero the smaller these should be. In non linear problems, however, since higher order terms exist, 
the accuracy need not improve indefinitely in the vicinity of a bifurcation point. 

As a continuation process starts and as long as the convergence rate is satisfactory, one should 
use <p = tii — «o in the global step. When approaching a singularity of the augmented system, 
for example, a bifurcation point (perfect or imperfect), it is necessary to approximate the nearly 
singular eigenfunctions and use them in the global step. An algorithm for approximating the almost 
singular subspace is therefore required for the method described here to be effective. It is presented 
in sec 6.2.1. 

The more serious difficulty which is encountered near singular points is related to badly ap- 
proximated smooth components on the coarse levels. Let <p h be a smooth eigenfunction of a linear 
operator L h . Because of smoothness it is also an approximate eigenfunction for the operator L ih 
which approximates L h . It can easily be shown [BT] that the rate of convergence of a two grid 
cycle using the operators L h and L 2h is approximately given by the quantity 


1 - 


X h 

A 2 * 


( 6 . 1 ) 


where X h ,\ ik are the eigenvalues of L h ,L 2h . When X h is close to zero, the above quantity may 
be close to or larger than one, i.e., the multigrid algorithm will have slow convergence or even 
divergence. Note the the number of nearly singular eigenfunctions is small ( compared to the total 
number of eigenfunctions) . For the bifurcation problems discussed in this paper, the number is 
small (just a few nearly singular components are possible). 

In nonlinear problems when the Frechet derivative of the operator is nearly singular, a similar 
difficulty occurs. An algorithm that is well suited for continuation near singular points is described 


next. 
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6.2 Modified MultiGrid Method 

As before consider a sequence of grids < M ) with mesh sizes h jt satisfying = h jt . 

Suppose on each grid an operator L* is given in such a way that L k (k < M) is an approximation 
to L* +1 . Assume that corrected coarse grid equations are used for levels k < l and that the 
subspace of badly approximated components, Ho , which is spanned by the set . . . <f> n } has 
been approximated (see section 6.2.1). The coarse grid equations are obtained in an analogue way 
the one obtained in [BT] for nearly singular problems. For k < l the equations to be solved for 


u*, (j = 1, . . . n) are given by 

L k (u k > X) = F k + f2n^i ( 6 . 2 ) 

i= i 

n 

N k (u k ,\, s ) = G k + J2 *?}■ wj (6.3) 

i - 1 

< «*> <t>) >= <*) + r) k a ) 1 <j<n, (6.4) 

where 

# = ^ _1 {r£ + i(u k+1 + X) - r k +l (n k+ \X)} + l£ +1 i» k+l ( k < l ) (6.5) 

"* = r*+l(# +1 ,A) (*</) (6.6) 

«- +1 = ^ +1 = 0 (j = !»•••«) (6-7) 

r‘+i(«* +1 ,A) = L k (7 k k+l u k+1 ,X) - I k +1 L k+1 (u k+ \ X) (6.8) 

^ + i(« fc+1 ,A) = N k (7 k +l u k+ \X,s) - N k+1 (u k+l ,X,s) (6.9) 

^ = n + i4>1 +1 U = 1 ,...«) ( 6 . 10 ) 

= L k (I k +1 u k+1 ,X) + I k +1 (F k+i - L k+1 {u k+1 , A)) (6.11) 

G* = W*(7* +1 u* +1 , A, s) + (G* +1 - JV fc+1 (u* +1 , A,s)) (6.12) 

= F* + £ (it < /) (6.13) 

i=i 

G* = G* + ^ (k < /) (6.14) 

p l + l — 


(6.15) 
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°) =< Ik+iu k+1 ,<t>j 


> + 


0 

a * +1 - < u* +1 ,^} +1 > 


if fc = / 
otherwise 


& i =ff i + ^ a 5 U = 1 »--- n ) 

a) =< <t> k j,4>) > (j = l,...n) 


(6.16) 

t&AT) 

(6.18) 


where r) k ,u k are the current approximations to fj k , U k respectively. Initial approximations are 
r) = 0, u k = l£ + iu* +1 . The interpolation in this case takes the form 


u k «— u k + j [tx* 1 — V] (6.19) 

+ (6.20) 

£*<-** + a V“ x (6.21) 

n 

F k + ^2 (6.22) 

JM-1 

G fc <- G* + ^2 ( 6 - 23 ) 

y=i 


6.2.1 Algorithm for Approximating the Nearly Singular Subspace 

As has been stated before, one needs approximations to the nearly singular components since these 
are slow to converge in the local relaxation process and are needed also in the formulation of the 
coarse grid equation . Actually what we compute are the nearly singular eigenfunctions of L u . 
The algorithm for doing that is an FMG algorithm. It begins by identifying these functions on 
the coarsest level. The device for that is Kaczmarz relaxation. It has the property of damping the 
error components which belong to eigenfunctions with relatively larger (in magnitude) eigenvalues. 
The ones that correspond to nearly singular components are the slowest to converge. By starting 
with a random guess W and relaxing the equation 


L u {u,\)W = 0 (6.24) 

enough times until the convergence rate becomes very slow, the approximation W, to the actual 
solution W = 0, lies mainly in the subspace of nearly singular functions. Starting with this as an 
approximation we can solve the following eigenvalue problem 
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L u {u,\)W - txW = 0 (6.25) 

H^ll 2 = 1 (6.26) 

in order to define a function in the almost singular subspace. Assume for the moment that this 
subspace is of dimension 1, and we solve for that subspace on the coarsest grid. A relaxation for 
that purpose consists of local and global steps. The local one is a Kaczmarz relaxation. The global 
one is of the form 

W +-WP (6.27) 

<L u (u,X)W,W> 

<W,W> 

where /? is such that the norm requirement is satisfied. When more that one grid is involved the 
relaxation is only the local one for the fine grids, while on the coarsest one also the global step is 
performed. Coarsening is done using FAS for the whole system (6.25)-6.26). The algorithm used 
is the FMG algorithm. 

In the event that the nearly singular subspace is of a higher dimension than one, that algorithm 
will converge slowly very soon; in fact, it can be observed even on the coarsest level when solving 
the discretized equations for (6.25)-(6.26). In case on any level (k > 1) a slowness of convergence 
is encountered, new eigenfunctions have to be introduced into the process, where at each time the 
initial approximation for such an eigenfunction is taken from the residuals, of the problems already 
in process, which have been orthogonalized to previous eigenfunctions found already, and for which 
an FMG algorithm is applied. The process should include a Ritz projection as in [BMR] for that set 
of functions. In this way an orthogonal set of functions that spans the nearly singular subspace is 
found. Each of these solves a system like (6.25)-(6.26) with a possibly different jx. When the finest 
grid is reached with the FMG algorithm, all the nearly singular eigenfunctions have been found. As 
mentioned before the accuracy needed for these functions depends on the closeness of the operator 
L u to singularity. A feasible approach is then to improve the accuracy of these eigenfunctions by 
one cycle of the eigenvalue problem per one cycle of the original problem in case convergence of 
the original problem remains slow. 
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6.2.2 Relaxation for Modified Equations 

As before when a system of equations has local and global unknowns, the relaxation is done in two 
steps. Here the local step is done for equation (6.2) keeping A and the right hand side fixed. This 
is done on all levels while the global step that is described next is performed for levels k <l g only, 
where l g is usually 1 or 2. Given {<£i, . . . , </>„} and {t/q, . . . , t/» n ) the global step is 


u *«- «* + £&*} (6.29) 

J=1 

Vj *- flj + Vj, 0 <j<n (6.30) 

\*-\ + 6 (6.31) 

where the unknowns {5, fa , . . . , /?„, »)i, . . . , r)„} are a solution to the following set of equations 

< L*(u* + £ Pi*)> X + «)» h >=< ~ Fk + £ >, 1 < J < n (6.32) 

i=i ?=i 

AT fc (u fc + f; A + 5, s) = G* + ^ f>juj (6-33) 

i=i y=i 

n 

< + II Pit), <t> k j >= + Vi ak 1 < 3 < n. (6.34) 

j=i 


The first of these equations can be approximated by using a Taylor expansion up to quadratic 
terms. Note that the linear terms only are not enough since they vanish at a singular point. 

Modified Multigrid Cycle (MMG) 

Having a relaxation procedure we can now formulate a modified multigrid method. For level 
k = M it is identical to MG; therefore, we describe it for levels k < l. It is denoted by 

(1 fc,u*4*,F*,G*,£*) <- MMG{k,u k ,^ k ,F k ,G k ,a k ) k<l (6.35) 

and is defined recursively as follows: 

IF k=l THEN solve (6.2)-(6.4) by enough relaxations (to achieve a desired accuracy) 

ELSE 

• Perform iq relaxation sweeps of section (6.2.2) on (6.2)-(6.4) resulting in a new approximation u k . 

• Starting with u* -1 = 7* -1 u* q* -1 = 0 ,F k ~ 1 ,G k ~ 1 ,^ c ~ 1 given by (6.13)-(6.17), make 
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7 successive cycles of the type 

4 - MMG(u k -\^ k - 1 ,F k - 1 ,G k - 1 ff k - 1 ). 

• Perform the corrections u k = u k + /jjL-Jtt* -1 — ,rf = rf + ,£* = £* + a k rf~ 1 

F k = F k + E"=i tjJ- VJ , G k = G k + Ei=i Vj ~ l «}. 

• Perform 17 relaxation sweeps on (6.2)-(6.4) starting with u* and yielding the 
final result of (6.35) 


The cycle with 7 = 1 is called V cycle or V(i7, 1 / 2 ), and the one with 7 = 2 is called W cycle or 
W{y u v 2 ). 

6.3 Modified Full MultiGrid Continuation Method (MFMGC) 

Given a sequence of problems 

L k {u k ,X) = F k (6.36) 

N k (u k ,X,s) = G k (6.37) 

that are defined on a sequence of grids Cl k with mesh sizes hie(k < M) and such that coarse grid 
operators are approximations to finer ones as before. We describe now a modified continuation 
multigrid method (MFMGC). It is assumed that two previous solutions (tii,Ai) and (tio,Ao) are 
given. 

continuation algorithm MFMGC 

• Choose As (for example «i — so ) and set s = s\ + As. 

Approximate the nearly singular subspace by 1-FMG. 

FOR k — 1, . .., M DO: 

• Perform the cycle ( u k ,ri k ,F k ,G k ,o r*) <— MMG(u k ,rf , F k ,G k a_ k ) N times. 

• If k < M calculate u fc+1 <— u fc+1 + Il£ +1 u fc . 

Call the resulting approximation on level M (u M , X) . 

• FORk=l, DO: 

u k = I* +1 u fc+1 
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min /x 2 subject to (6.38) 

{ 0 > s ) 

< 2/(u + /3<f> , A + £),<£ >=< jFi, <f> > (6.39) 

< Z/ U (tx 4* A + 5)^> — fi < <f>,(f> >=< F 2 ><f> > . (6.40) 

This by itself is not enough since the second constraint equation (3.12) will converge very slowly in 
the <f> direction. This can be taken care of by changing /x on the coarsest levels such that equation 
(6.40) holds there. Note that a change of /x in this way is effectively relaxing this equation in 
the <f> direction. This step on the coarsest levels is preceded by updating the norm of <f> according 
to its constraint by multiplying it by the appropriate constant. To summarize the algorithm is 
basically an FMG algorithm using the FAS formulation where the local relaxation is combined 
with two global steps. One is performed on the finest level and the other on the coarsest levels. It 
is demonstrated in section 7 that 2-FMG-W cycles solve for the location of a bifurcation point, to 
the level of discretization error. The minimization problem (6.38)-(6.39) is solved by approximating 
it by a minimization problem obtained by expanding L(u + /?<£, A + 6) and L u (u + A + 6) <f> by 
Taylor series taking up to quadratic terms in the first one and only linear terms in the second; that 
is, 


min/x 2 subject to (6.41) 

( 0 » 5 ) 

/9A + 5B + lp 2 C + p6D + U 2 E=< F t - L(u, \),<f> > (6.42) 

z z 

PC + SD — fi < <f>, (f> >=< Fi — L u (u, X)<f>,(f>> (6.43) 

where 


A =< L u (u,X)4>,<f> >, B =< L\(u,X),<f> > C =< L u u(u,X)(f><f>,<f> > 


D =< <}> >, E=<L\\,<f>>. 
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7 Numerical Experiments 

In this section we demonstrate some of the algorithms presented in the paper. We consider partial 
differential equations of the form 

Ati + /(x,y,u,A) = 0; (x, y) € (0, 1) x (0, 1). (7.1) 

The discretization used for all cases is the standard five point Laplacian for the A u term in the 
equation, and a pointwise approximation to f(x,y,u, A). The grids used were equally spaced in 
both directions. 

7.1 Locating Limit Points 

In order to demonstrate the algorithm for locating limit points we, have chosen the well known 
Bratu problem. That is, 

/(x,y,u,A) = Ae u . (7.2) 

A bifurcation diagram for this problem is given schematically in Figure 1. 



Table 1 shows the results of a 2-FMG W cycle. The coarsest grid has a mesh size of ho = .25. 
Bilinear interpolation used for 7*_ x , bicubic for II*_ X , and a 9-point full weighting operator for 



29 


l£ x . The equations solved on each level are of the form 


L k (u k ,X k ) = F k 

(7.3) 

L k {u k ,X k )<p k = G k 

(7.4) 

IMI 2 = r* 

(7.5) 


where F M = 0 ,G M = 0, r m = 1 and for k < M,and F k ,G k ,r k for k < M are obtained from the 
FAS equations. The relaxation used was Newton Gauss- Seidel for u k in the first equation , Gauss 
Seidel for <p k in the second ( keeping u k ,X k fixed ). This local relaxation was employed on all levels. 
On the coarsest level also a global step was used as described in section 5.3. The residuals given 
in the table are for the currently finest level at the end of the cycle. The two residuals at each row 
(starting at the left one) are for equation (7.3), (7.4) respectively. In the notation of section (5.1) 
we have used v\ = 2 and U 2 = 1. It is clearly seen that 1-FMG solves for the location of the limit 
point to the level of discretization errors. The approximation to the limit point X k obtained by 
1-FMG converges as 0(h 2 ). The same rate is obtained for < L k , <p h >. 

7.2 Continuation near Singular Points 

Continuation in the vicinity of singular points should be done with a special multigrid algorithm 
like MFMGC. We show here the difficulties encountered with a standard method and compare the 
results to the modified algorithm MFMGC. Tables 2 and 3 show a two grid cycling in the vicinity 
of a nearly singular point of the differential equation. In this case the differential equation is given 
by (7.1) with /(x,y, u,A) given by 

/(x,y,u,A) = A(u + u 2 ) + e. (7.6) 

The choices e = 0 and e = —.01 are shown schematically in Figure 2. 

The numerical experiments in this section were for the case e = - .01. The reason for using only 
a two grid experiment was to see clearly the fact that the FMGC algorithm has difficulties while 
MFMGC has a very good convergence rate. Note that even in case the coarsest grid problem is 
solved exactly the two grid cycle has a very slow convergence rate which indicates that the coarse 
grid does not approximate well enough some components. The experiment is given for a value of A 
which is in a vicinity for which the linearized problem is almost singular. 
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Figure 2: Bifurcation diagram 

7.3 Locating Bifurcation Points 

In this section the algorithm presented in section (6.4) is demonstrated. Two cases are given 
for which the differential equation has bifurcations at A = 7r 2 (m 2 + n 2 ) where n,m are integers. 
A simple way of constructing a problem in the form of (7.1) for which bifurcations exist at the 
above mentioned A’s [K] is to begin by prescribing a branch of solutions of the form ti(x,y,A) = 
g(A)C/o(£jy) where y(A), Uo(x,y) are given. This is achieved if 

f{x,y,u(x,y, A), A) = Au(x,y,X). 

In order that branches will bifurcate from this branch at the above A ’s we require that 

fu[x,y,u{x, y,A),A) = 1. 

As a special case of this we take 

f(x,y,u,X) = -q(X)AU 0 (x,y) + XP(u - q(\)U 0 (x,y))] P(z) = z + z 2 . (7.7) 


We consider two cases: (i) Uo(x ) y) = x(x — l)y(y- 1), and (ii) J7o(x,y) = sm(7rx)sin(7ry). The first 
is a case of a perfect bifurcation for the dicrete levels while the other is an imperfect bifurcation 
for the discrete levels. Tables 3 and 4 show the result of the 2-FMG-W algorithm for locating 
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the first bifurcation point. Intergrid transfers used are identical to the ones used in section 7.1. 
The two residuals (left and right) are for the first and second constraint equations (3.11) ,(3.12) 
respectively. The approximated singular point Aj converges at a rate which looks as 0(h) for the 
imperfect bifurcation case and as 0(h 2 ) for the other case. Also note that < L%,<p h > in case 
(ii) converges to zero as O(h). As the actual value of this quantity shows, it is not enough to 
consider single grid experiments when locating bifurcation points. A refinement is essential for 
distinguishing between actual limit points and imperfect discrete bifurcation. The results of both 
tables clearly demonstrate the effectiveness of the algorithm described. 
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level 

cycle # 

|| residuals ||2 

A 

< > 

i 

5 

0.332e-14 

0.435e-14 

6.69051 

1.86767 

2 

1 

0.390e-l 

0.111e+l 

6.79469 

1.98973 


2 

0.172e-2 

0.123e-l 

6.78333 

1.99357 

3 

1 

0.135e-l 

0.417e+0 

6.80282 

2.02417 


2 

0.308e-3 

0.136e-l 

6.80217 | 

2.02440 

4 

1 

0.42 le-2 

0.130e+0 

6.80669 

2.03212 


2 

0.129e-3 

0.213e-2 

6.80665 

2.03214 

5 

1 

0.180e-2 

0.351e-l 

6.80776 

2.03408 


2 

0.603e-4 

0.209e-3 

6.80775 

2.03408 


Table 1: Locating limit point 



|| residuals H 2 

A 

0.833e-16 

18.99997 

0.165e-01 


0.108e-01 

0.101e-01 

18.99997 

0.956e-02 


0.517e-02 

0.804e-02 

19.00000 

0.792e-02 


0.310e-16 

0.142e-01 

19.00002 

0.141e-01 


0.393e-16 

0.117e-01 

19.00003 

0.117e-01 


0.555e-16 

0.721e-02 

19.00004 


|| residuals H 2 

A 

0.500e-16 

18.99999 

0.374e-02 


0.500e-16 

0.795e-03 

18.99999 

0.342e-03 


0.340e-16 

0.544e-04 

18.99999 

0.187e-04 


0.139e-16 

0.296e-05 

18.99999 

0.132e-05 


0.572e-16 

0.178e-06 

18.99999 

0.116e-06 


0.310e-16 

0.140e-07 

18.99999 


Table 2: Algorithm FMGC 


Table 3: Algorithm MFMGC 
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level 

cycle # 

|| residuals ||a 

A 

< L x ,4>> 

m 

5 

0.915e-10 

0.176e-7 

18.74516 

-6.58175e-8 

2 

1 

0.264e-2 

0.381e+0 

19.28335 

3.64700e-3 


2 

0.725e-3 

0.387e-l 

19.46164 

4.43231e-4 

3 

1 

0.509e-3 

0.264e-l 

19.67002 

1.01890e-4 


2 

0.364e-4 

0.336e-l 

19.67527 

1.04083e-5 

■ 

1 

0.715e-4 

0.280e-2 

19.72318 

2.98171e-6 

■ 

2 

0.235e-5 

0.110e-3 

19.72336 

-3.31378e-8 

5 

1 

0.133e-4 

0.297e-3 

19.73545 

-3.63521e-6 


2 

0.97 4e-7 

0.129e-3 

19.73528 

-7.34759e-7 


Table 4: Locating a bifurcation point: Perfect bifurcation 


level 

cycle # 

|| residuals ||2 

A 


m 

5 

0.186e-13 

0.628e-13 

9.48662 

0.56972 

2 

1 

0.160e-l 

0.545e+0 

13.46777 

0.20472 


2 

0.291e-4 

0.125e-l 

13.46833 

0.20487 

3 

1 

0.308e-l 

0.901 e+0 

16.21966 

0.85237e-l 


2 

0.149e-2 

0.400e-l 

16.21687 

0.85390e-l 

■ 

1 

0.510e-2 

0.591e+0 

17.86465 

0.38793 e-1 

■ 

2 

0.113e-3 

0.135e-l 

17.86442 

0.38813e-l 

5 

1 

0.829e-3 

0.335e+0 

18.77085 

0.18483e-l 


2 

0.542e-5 

0.421e-2 

18.77096 

0.18484e-l 


Table 5: Locating a bifurcation point: Imperfect bifurcation 
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